{"nbformat":4,"nbformat_minor":0,"metadata":{"colab":{"provenance":[],"collapsed_sections":["p0jrcMSbKdix","CVwOxESDKpS_"],"authorship_tag":"ABX9TyMQp5RUAOxIUJPBieHCsO9/"},"kernelspec":{"name":"python3","display_name":"Python 3"},"language_info":{"name":"python"}},"cells":[{"cell_type":"code","source":["# Install packages\n","import pandas as pd\n","import numpy as np"],"metadata":{"id":"YqgTZcJn3bP2"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Mount my Google Drive\n","from google.colab import drive\n","drive.mount('/content/drive')\n","\n","\n","# NOTE paths below are specific to my computer/Google Drive and will need to be\n","# corrected for the code to run."],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"kXDXoLGtKEk9","executionInfo":{"status":"ok","timestamp":1715628856861,"user_tz":420,"elapsed":1210,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"d523f769-75da-4e74-82da-913cac546e77"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["Drive already mounted at /content/drive; to attempt to forcibly remount, call drive.mount(\"/content/drive\", force_remount=True).\n"]}]},{"cell_type":"markdown","source":["# Electric Load"],"metadata":{"id":"p0jrcMSbKdix"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","elec_use_daily_person = 1.5 # (see lit notes)\n","household_size = 8          # 8.3 average for Senegal (see lit notes)"],"metadata":{"id":"PwzUTDOAH9_2"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Path to load profile in Google Drive\n","load_path = '/content/drive/My Drive/1 - projects/Senegal USAID/data/Load/Huld 2017 Fig. 4 Low nighttime consumption profile.csv'"],"metadata":{"id":"6JjW-xZL15Ez"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Import load data from csv of extracted profile\n","data = np.genfromtxt(load_path, delimiter=',')\n","\n","x = data[:, 0]  # First column (index 0) is x (hr)\n","y = data[:, 1]  # Second column (index 1) is y (kW)\n","\n","# Check that they contain the same number of points\n","len(x) == len(y)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"RlqYv6OgN7JA","executionInfo":{"status":"ok","timestamp":1715628856861,"user_tz":420,"elapsed":7,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"6f75f8d2-ae97-4453-e17e-fe20100b7db1"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["True"]},"metadata":{},"execution_count":83}]},{"cell_type":"code","source":["# Convert load data to a dataframe\n","df = pd.DataFrame(data, columns=['x', 'y'])\n","\n","# Generate integer x values from 1 to 24\n","# and convert to float for consistency with df column of x values\n","all_x_values = pd.Series(range(1, 25), name='x').astype(float)\n","\n","# Merge the given x values with all_x_values\n","merged_df = pd.merge(df, all_x_values.to_frame(), on='x', how='outer')\n","\n","# Sort to get x-values in sequential order\n","merged_df = merged_df.sort_values('x')\n","\n","# Perform linear interpolation to fill NaN values for y\n","merged_df['y'] = merged_df['y'].interpolate(method='linear')\n","\n","# Extract the power consumption number, y, for integer x values from 1 to 24.\n","merged_df = merged_df[merged_df['x'].isin(all_x_values)]\n","\n","# Move 24 (midnight) value to beginning of time series\n","# and re-index\n","# and set the x-value to 0\n","midnight_row = merged_df.iloc[-1:]\n","other_rows = merged_df.iloc[:-1]\n","loadprofile_df = pd.concat([midnight_row, other_rows])\n","loadprofile_df = loadprofile_df.reset_index(drop=True)\n","loadprofile_df.at[0,'x'] = 0.0\n","\n","# Convert column x to integers\n","loadprofile_df['x'] = loadprofile_df['x'].astype(int)"],"metadata":{"id":"iTDpqhWPymik"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Verify that it is scaled to 300 kWh/day consumption per Huld (2017).\n","original_scale = sum(loadprofile_df['y'])\n","\n","original_scale"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"oDOjptFZGid3","executionInfo":{"status":"ok","timestamp":1715628856861,"user_tz":420,"elapsed":6,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"1f967880-598e-457b-9b45-c329c473dc26"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["299.7066244563945"]},"metadata":{},"execution_count":85}]},{"cell_type":"code","source":["# Rescale to Senegal household consumption\n","scaling_factor_household = elec_use_daily_person * household_size / original_scale\n","loadprofile_df['y'] *= scaling_factor_household"],"metadata":{"id":"X7wTHh0FIMvb"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Verify that it is scaled for household to elec_use_daily_person * household_size.\n","new_scale = sum(loadprofile_df['y'])\n","\n","print(new_scale)\n","print(elec_use_daily_person * household_size)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"rAcvZQMYJMEN","executionInfo":{"status":"ok","timestamp":1715628856861,"user_tz":420,"elapsed":5,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"579f1ad8-af44-420d-a7eb-e35951c4f322"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["11.999999999999998\n","12.0\n"]}]},{"cell_type":"code","source":["# Make into a yearly load profile by repeating 365x.\n","yearly_loadprofile = np.tile(loadprofile_df.values, (365, 1))\n","yearly_loadprofile_df = pd.DataFrame(yearly_loadprofile)\n","\n","# Name columns\n","# assuming power value is constant through the hour, so that demand = kW * 1h\n","yearly_loadprofile_df.columns = ['Hr of day','Demand (kWh)']"],"metadata":{"id":"ZXSgClae__Mt"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# PV System"],"metadata":{"id":"CVwOxESDKpS_"}},{"cell_type":"code","source":["# VARIABLES\n","\n","size_kWp = 3.5"],"metadata":{"id":"pgM6yBHyWa6v"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","sim_size_kWp = 4        # PVWatts default\n","panel_eff = 0.19    # PVWatts default for Crystalline Silicon\n","tilt = 15           # close to latitude\n","orientation = 180   # S"],"metadata":{"id":"0tgSOe41K166"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Calculate array area\n","PV_area = size_kWp / panel_eff # kW * (1 m2/kW under test conditions) * kW/kW\n","PV_area # (m2)"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"fDUKA-17LgtW","executionInfo":{"status":"ok","timestamp":1715628856862,"user_tz":420,"elapsed":5,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"23de847d-018e-48bd-fafc-89157bcdf996"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["18.42105263157895"]},"metadata":{},"execution_count":91}]},{"cell_type":"code","source":["# Path to PV output profile in Google Drive\n","PVdata_path = '/content/drive/My Drive/1 - projects/Senegal USAID/data/PV/pvwatts_hourly_Dakar_4kW_rooftop.csv'"],"metadata":{"id":"IMdUnOS7BIh6"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Import load data from csv of extracted profile\n","PVoutput_df = pd.read_csv(PVdata_path, skiprows=31)"],"metadata":{"id":"toKTNbQSBIkM"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Scale output according to PV system size\n","PVoutput_df['AC System Output (W)'] = PVoutput_df['AC System Output (W)'] * size_kWp / sim_size_kWp\n","PVoutput_df['DC Array Output (W)'] = PVoutput_df['DC Array Output (W)'] * size_kWp / sim_size_kWp"],"metadata":{"id":"Dj-C9yPZW74H"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Convert AC system output (W) values to kWh values for the previous hour\n","# by converting W to kW and assuming power provided is constant\n","# through the previous hour, so that PV output is in kW * 1h\n","\n","PVoutput_df['PV Output (kWh)'] = PVoutput_df['AC System Output (W)'] / 1000"],"metadata":{"id":"PpUKLStTDvVR"},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":["# Combine working data into 1 dataframe and check reliability"],"metadata":{"id":"10wChh4fDAeN"}},{"cell_type":"code","source":["# Filter load dataframe\n","yearly_loadprofile_df_filtered = yearly_loadprofile_df[['Demand (kWh)']]"],"metadata":{"id":"Kh1q56W_DG8k"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Filter PV output dataframe\n","PVoutput_df_filtered = PVoutput_df[['Month','Day','Hour','PV Output (kWh)']]"],"metadata":{"id":"jdT7U50wDpeK"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Concatenate into 1 df\n","df = pd.concat([PVoutput_df_filtered, yearly_loadprofile_df_filtered], axis=1)"],"metadata":{"id":"vLP-ENh8EEEE"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Compare total PV output to total demand.\n","PVoutput_to_demand_annual = sum(df['PV Output (kWh)'])/sum(df['Demand (kWh)'])\n","PVoutput_to_demand_annual"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"dCOBwUBALyIs","executionInfo":{"status":"ok","timestamp":1715628857095,"user_tz":420,"elapsed":4,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"29ba60cc-9ea0-4d65-e839-29ce13479160"},"execution_count":null,"outputs":[{"output_type":"stream","name":"stdout","text":["5585.278909874995\n","4379.999999999882\n"]},{"output_type":"execute_result","data":{"text/plain":["1.2751778333048278"]},"metadata":{},"execution_count":99}]},{"cell_type":"markdown","source":["# Battery and check reliability"],"metadata":{"id":"aPs432oDKpyM"}},{"cell_type":"code","source":["daily_demand = sum(df['Demand (kWh)'])/365\n","daily_demand"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"1lGA_mNkaCLN","executionInfo":{"status":"ok","timestamp":1715628857095,"user_tz":420,"elapsed":3,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"20e7a88c-b59f-47b3-bfa5-388f44aa7245"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["11.999999999999677"]},"metadata":{},"execution_count":100}]},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","capacity = 6.5          # (kWh)\n","\n","efficiency = 0.95     # round-trip efficiency\n","\n","charge_rate = 0.8     # (capacity/hr)\n","discharge_rate = 1.0  # (capacity/hr)\n","\n","reliability = 0.95    # portion of hours that demand should be met by PV + batt"],"metadata":{"id":"Ymqcg8YfTT8S"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["class Battery:\n","    def __init__(self, capacity, efficiency, init_charge=0):\n","        \"\"\"Initialize battery parameters.\"\"\"\n","        self.capacity = capacity    # Total energy storage capacity (kWh)\n","        self.efficiency = efficiency  # Round-trip efficiency (0.0 - 1.0)\n","        self.charge = init_charge     # Initial state of charge (kWh)\n","        self.max_charge_rate = 0.5 * capacity  # Maximum charge rate (kW)\n","        self.max_discharge_rate = 0.5 * capacity  # Maximum discharge rate (kW)\n","\n","    def update_charge(self, demand, pv_output):\n","        \"\"\"\n","        Update the battery charge based on demand and PV output.\n","\n","        Args:\n","            demand: Current electricity demand (kW)\n","            pv_output: Current PV power output (kW)\n","        \"\"\"\n","\n","        if net_energy > 0:  # Surplus energy\n","            max_charge = min(net_energy, self.max_charge_rate, self.capacity - self.charge)\n","            self.charge += max_charge * self.efficiency\n","        else:  # Energy deficit\n","            max_discharge = min(-net_energy, self.max_discharge_rate, self.charge)\n","            self.charge -= max_discharge / self.efficiency\n","\n","        # Ensure charge stays within bounds\n","        self.charge = max(0, min(self.charge, self.capacity))"],"metadata":{"id":"zfzUS3ceQkHv"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Instantiate battery\n","battery = Battery(capacity=capacity, efficiency=efficiency)"],"metadata":{"id":"otfwk1p_amVI"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Simulate battery operation over each row (hourly timestep)\n","for index, row in df.iterrows():\n","    battery.update_charge(row['Demand (kWh)'], row['PV Output (kWh)'])\n","    df.at[index, 'Charge (kWh)'] = battery.charge"],"metadata":{"id":"F45pRnOjPWiH"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Calculate battery output\n","df['Previous Charge (kWh)'] = df['Charge (kWh)'].shift(1).fillna(0)\n","\n","df['Battery Output (kWh)'] = df['Previous Charge (kWh)'] - df['Charge (kWh)']"],"metadata":{"id":"QyDhyqPSSKCE"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Identify hours where PV + battery combined do not meet the demand\n","df['Insufficiency'] = df['Demand (kWh)'] > (df['PV Output (kWh)'] + df['Battery Output (kWh)'])\n","\n","# Quantify unmet demand in each hour\n","df['Unmet Demand (kWh)'] = df['Demand (kWh)'] - df['PV Output (kWh)'] - df['Battery Output (kWh)']\n","df['Unmet Demand (kWh)'] = df['Unmet Demand (kWh)'].clip(lower=0)\n","\n","# Check whether it meets reliability criteria\n","sum(df['Insufficiency']) < (8760 * (1-reliability))"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"oQCPRKEeXgdT","executionInfo":{"status":"ok","timestamp":1715628858209,"user_tz":420,"elapsed":7,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"b507c62c-f756-46ac-bd38-cf17f19cdc1e"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["True"]},"metadata":{},"execution_count":106}]},{"cell_type":"code","source":["# Quantify reliability\n","hours_met = (1-sum(df['Insufficiency']))/8760\n","(1 - sum(df['Insufficiency'])/8760)*100"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"00IZXe-r6hRk","executionInfo":{"status":"ok","timestamp":1715628858209,"user_tz":420,"elapsed":6,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"c4bf1b66-6c63-4c78-d09d-b653add3c97a"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["95.52511415525115"]},"metadata":{},"execution_count":107}]},{"cell_type":"markdown","source":["# Cost Estimate"],"metadata":{"id":"SpqgrioGKqFq"}},{"cell_type":"code","source":["# ASSUMPTIONS\n","\n","PV_cost_per_kWp = (950 + 1150)        # USD, from Szabo et al. 2021, including modules + installation and BOS\n","battery_cost_per_kWh = 400            # USD, from Szabo et al. 2021, for Li-ion\n","\n","# NOTE\n","# O&M costs (for battery and panel maintenance) are ignored here!"],"metadata":{"id":"NYBwZpfJ3N4r"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Estimate PV system cost\n","PV_cost = PV_cost_per_kWp * size_kWp\n","PV_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"9yR0xWBD5zOC","executionInfo":{"status":"ok","timestamp":1715628858209,"user_tz":420,"elapsed":5,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"4375e970-1afc-45fd-92cd-e39327874d54"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["7350.0"]},"metadata":{},"execution_count":109}]},{"cell_type":"code","source":["# Estimate battery cost\n","battery_cost = battery_cost_per_kWh * capacity\n","battery_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"HKq1Bby06H51","executionInfo":{"status":"ok","timestamp":1715628858209,"user_tz":420,"elapsed":4,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"f7a8e2c1-82ac-4756-eb90-09723465c684"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["2600.0"]},"metadata":{},"execution_count":110}]},{"cell_type":"code","source":["# Estimate total system cost\n","total_cost = PV_cost + battery_cost\n","total_cost"],"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"id":"0IG0uJ0b6ICf","executionInfo":{"status":"ok","timestamp":1715628858209,"user_tz":420,"elapsed":4,"user":{"displayName":"Amanda Smith","userId":"00495916758981390713"}},"outputId":"479a2a85-6a93-45fd-8335-71531120a40e"},"execution_count":null,"outputs":[{"output_type":"execute_result","data":{"text/plain":["9950.0"]},"metadata":{},"execution_count":111}]},{"cell_type":"markdown","source":["# Output file"],"metadata":{"id":"O-Of6An97r1S"}},{"cell_type":"code","source":["# Create output df\n","output_data = {\n","    'People': [household_size],\n","    'PV capacity (kWp)': [size_kWp],\n","    'PV panel area (m2)': [round(PV_area, 2)],\n","    'PV system cost ($)': [PV_cost],\n","    'Battery capacity (kWh)': [capacity],\n","    'Battery cost ($)': [battery_cost],\n","    'Total cost ($)': [total_cost]\n","}\n","output_df = pd.DataFrame(data=output_data)"],"metadata":{"id":"1wYfFM0o74LQ"},"execution_count":null,"outputs":[]},{"cell_type":"code","source":["# Create output file\n","output_file_name = f\"results_tilt{tilt}deg_orientation{orientation}.csv\"\n","output_file_path = \"/content/drive/My Drive/1 - projects/Senegal USAID/\"\n","output_df.to_csv(output_file_path + output_file_name, index=False)  # exclude row numbers"],"metadata":{"id":"nNQDw7zoSr2t"},"execution_count":null,"outputs":[]}]}